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Abstract We study numerical simulations of large {N « 
10^) two-dimensional quasi-static granular assemblies sub- 
jected to a slowly increasing deviator stress. We report 
some peculiarities in the behavior of these packings that 
have not yet been adressed. The number of sliding con- 
tacts is not necessarily related to stability: first the num- 
ber of sliding contacts rises linearly and smoothly with 
the applied stress. Then, at approximately half the peak 
stress, the increase slows down, a plateau develops, and a 
decrease follows. The spatial organization of sliding con- 
tacts also changes: during the first half of the simulation, 
sliding contacts are uniformly distributed throughout the 
packing, but in the second half, they become concentrated 
in certain regions. This suggests that the loss of homo- 
geneity occurs well before the appearance of shear bands. 
During the second half events appear where the number of 
sliding contacts drops suddenly, and then rapidly recovers. 
We show that these events are in fact local instabilities in 
the packing. These events become more frequent as fail- 
ure is approached. For these two reasons, we call these 
events precursors, since they are similar to the precursors 
recently observed in both numerical [H [2] and experimen- 
tal [21 m O [5] studies of avalanches. 



> 1 Introduction 



The principal goal of many numerical studies of quasi- 
static granular materials is to establish a connection be- 
tween the macroscopic, continuum-like response of the ma- 
terial, and its microscopic, grain- level state [7llHlin|- One 
hopes to explain the macroscopic behavior using a small 
number of microscopic parameters that could then be used 
in a constitutive model. The microscopic physical mean- 
ing of these parameters would provide the model with a 
physical justification, making it more general. Several can- 
didates have been examined, such as the fabric tensor TO] , 
force chains [UlIIl] or force distributions 13J, and sliding 
contacts [T3J[T5], In the limit of isostatic networks, stress 



paths have been calculated directly from the fabric tensor 
[16] . However, for arbitrary, frictional granular systems 
the subspace of all allowed force networks (and therefore 
stress paths) have to be included in the analysis [TTj . This 
leads to ambiguous values of the microscopic parameters. 
Unfortunately, the future evolution of the system depends 
on these parameters. 

Although the goal of a general, physically-based stress- 
strain relation has not been attained, several advances 
have been made in this direction. One important point 
is the role of sliding contacts. The contacts between the 
individual grains are assumed to be governed by Coulomb 
friction. Thus each contact can be in one of two states: 
either "sliding" or "non-sliding". The Coulomb friction 
law is sufficient to explain the incremental non-linearity of 
granular materials: previous studies [HI [H] have shown 
that increasing the load leads to an augmentation of the 
number of sliding contacts, which in turn causes the stiff- 
ness of the material to decrease. If the loading direction is 
reversed, the majority of the sliding contacts become non- 
sliding, leading to an abrupt stiffening. For this reason, it 
is believed that the contact status (sliding or non-sliding) 
is the most relevant microscopic variable. 

This connection between softening and sliding contacts 
leads to the expectation that the number of sliding con- 
tacts will continue to increase, right up to the time when 
the packing fails. Although this expectation seems rea- 
sonable, it has not been checked. Most numerical studies 
have focussed on cyclic loading far from failure [THl [20] > 
or the critical state [211 HI]. In this paper, we examine 
granular packings under increasing load up to the time 
of failure. We find that the expectation of a monotoni- 
cally increasing number of sliding contacts is false. The 
number of sliding contacts attains a clear maximum well 
before failure. Thus the density of sliding contacts cannot 
be used as an internal variable. Packings with the same 
number of sliding contacts may be in very different states. 

Our work also shows the existence of precursors in bi- 



axial tests. Previously, precursors of avalanches have 
been identified in both numerical [3 [5] and experimen- 
tal [51 mini [3] studies. In these studies, the inchnation of a 
static granular bed is slowly increased, until an avalanche 
occurs. Preceding the avalanche are numerous local re- 
organizations of the packing. These events become more 
and more frequent as the angle of inclination is increased. 

In biaxial tests, precursors also become more frequent 
when the failure is approached. They are triggered by 
changes in a localized region, and they produce sound 
waves that propagate outwards. At the origin of a pre- 
cursor there is always a local instability that triggers the 
precursor. These instabilities resemble those that trigger 
failure in very small packings |23j . The role of the pre- 
cursors in the failure of large assemblies has yet to be 
investigated. 

After a brief description of the simulation method, we 
introduce to the simulation parameters and define the 
quasi-static limit in Sec. 2] Next we give a description 
of the simulation in Sec. [3 We show how the stress-strain 



curve (Sec. 3.11 and the kinetic energy (Sec. 3.2) evolve 
as the system approaches the failure. Furthermore we ex- 
amine the volume, the injected power, and the number of 
contacts in Sec. 13.31 Thereafter we discuss in Sec. 13. 41 the 



evolution of the number of sliding contacts (Sec. 3.4.1 ), the 



average force transmitted at sliding contacts (Sec. 3.4.21, 
the contact status transitions (Sec. [3.4.3 and the spatial 
organization of sliding contacts (Sec. 3.4.4). Last but not 
least, we discuss the two regimes of qualitatively different 
granular behavior prior to the failure (Sec 3.5). In a fur- 
ther section (Sec. El, we examine a precursor carefully, an- 
alyzing the number of sliding contacts (Sees 4.2.1 4.2.4), 
the stress-strain relation (Sec. 4.3), the stability of the 



energy (Sec. 4.2.3) 



packing (Sec. 4.2.2) as well as the evolution of the kinetic 
After that we show that the sliding 

and 
We 



contacts tend to cluster when a precursor appears 
the clusters disappear again afterwards (Sec 



4.2.4) 



show that the precursor is localized in the packing, but 
that the vibrations that appear afterwards travel through 
the packing. The section is terminated by a summary of 



the precursor results in Sec. 4.4 We conclude our work in 
Sec. [5] by some speculations on the significance of precur- 
sors for the failure. 



2 Numerical procedure 

2.1 Contact model 

Grains are modeled as disks, and their interactions are 
calculated using the common "soft-sphere molecular dy- 
namics" method [M|. The force at the grain contact is 
generated by a linear dissipative spring whose length is 
given by the overlap distance £>„: 



-^n — '^n^n ^n-^n- 



(1) 



Here, kn is the length independent spring stiffness and 
the damping coefficient 7„ controls the energy dissipation. 
The overlap distance Dn is calculated from the radii r^ and 
rj of the touching particles and their positions x; and Xj , 



When the surfaces of the two touching disks move relative 
to each other a second force F^ arises, directed tangent 
to the particle surfaces. In analogy with the normal force 
defined in Eq. (fTj) we have 



Ft = -ktDt - jtDt 



(3) 



In our case kt = k„ and "ft = Jn- Determining the change 
in the tangential spring length Dt involves both transla- 
tional movement and rotation of the two touching particles 
i and j, 



Dt = -riUJi - rjUjj 



Ti + r 



r, - Dn 



V(v.-v,)-t. (4) 



^n \^i ^j I ^i ^j 



(2) 



Here, uji and ujj are the angular velocities of the touching 
particles and t is a vector tangent to the particle surfaces 
at the point of contact. The factor in front of the last term 
is needed to account for the overlap of the particles. 

In our simulation we allow only for repulsive forces, 
Fn > 0. We also enforce the Coulomb condition at each 
contact. 

fiF,, >\Ft\. (5) 

At a given moment, contacts where the strict inequal- 
ity holds are nonsliding contacts, whereas contacts with 
/iF„ = \Ft\ are sliding. When two particles separate, the 
contact "opens" , which can be considered as a third pos- 
sible contact status. In the following, we abbreviate the 
possible statuses of the contacts by O (open) , S (sliding) 
and C (closed or non-sliding). 

Contact status changes, i.e., transitions between "slid- 
ing" and "non-sliding" play an important role in this pa- 
per. A contact undergoes a transition C -^ S (non-sliding 
to sliding) when applying Eq. m would lead to a violation 
of Eq. ([5]). The reverse transition {S — > C) occurs when 
Dt in Eq. W changes sign or F„ increases. The transi- 
tions S ^ O and C ^ O occur when two touching grains 
separate, and the reverse transitions O -^ S" or O — >■ C 
occur when two grains come into contact. 

2.2 Units and parameters 

Three parameters are set to unity for all simulations: the 
particle density p, the initial system length L, and the 
pressure p. This defines our system of units. In two di- 
mensions, the unit of force is / = pL, and the unit of 
energy is E = pL^ . The unit of mass is m = pL"^ , whereas 
the time is measured in units of i = L^J p/p. The spring 
stiffness has a value of kn,t = 1600p. This leads to overlaps 
that are a small fraction of the radius: on average, we have 
F>n,t/FL « 0.3%. The small overlaps avoid the creation of 
many additional contacts. The damping is 7„_t = 0.19m/i, 
and the Coulomb friction coefficient is /i = 0.25. No grav- 
ity is applied to the particles. Unless otherwise mentioned, 
our systems have N = 16384 particles with average parti- 
cle mass m = 2.8 10~^/5L^. The total mass of the system 
is therefore of order unity. 

2.3 Boundary conditions 

We apply biaxial boundary conditions. These are easy 
to implement and simple to handle. In each direction. 
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Figure 1: Biaxial boundary conditions. The walls are 
smooth and can move normal to their surfaces. 



the granular packing is delimited by a light-weight (m — 
O.Olrfi) moveable wall parallel to one of the coordinate 
axes, as shown in Fig.jT] A force is applied to each bound- 
ary that can be constant or time-dependent. In this way, 
one can fix the average stress inside the granular pack- 
ing. The deformation of the packing can be determined 
by monitoring the movements of the boundaries. 

In our setup, the walls are perfectly slippery, i.e. they 
exert only normal forces. This has several advantages. 
First, the forces on opposite walls are guaranteed to be 
equal. Second, the number of degrees of freedom is re- 
duced. And last but not least, larger systems are more 
homogeneous. 

2.4 Preparation of initial conditions 

The 2 dimensional granular medium is made up of 
disks with radii r uniformly distributed within the range 
[O.Trniaxj Ij'max] with Tmax = 3.5 10^"^ and initial velocity 
in the range v^^Vy G [—0.5,0.5]. The initial radii and ve- 
locities are chosen using a random number generator, and 
a series of packings is generated by changing the "seed" . 
The resulting assemblies obtained essentially differ in con- 
tact topology and the spatial distribution of the grain 
sizes. In this paper, we study in detail the behavior of 
one packing containing 16384 disks. However, we com- 
pared the behavior of this packing to three other packings 
of this size, and we found that the results are always the 
same. Therefore the results do not depend on the config- 
uration but apply for all packings of this size. 

The frictionless particles are initially separated, but a 
constant external force /o = 0.8/ is applied to each of 
the walls, so that they move inward and compress the 
grains into a packing of approximate size 0.8L capable of 
transmitting normal forces. 

During the compression process, kinetic energy is re- 
moved by the damping at particle contacts. Certain mo- 
tions, however, require special care. For example, the ve- 
locity of the center of mass cannot be damped by con- 
tact forces, and thus a global viscous damping is applied 
for 10£. The kinetic energy in the assembly decreases 
to 2.6 10~^^E. The main reservoir of remaining kinetic 
energy are particles without contacts. To remove their 
energy, a viscous damping force opposing the individual 



grain movement is then applied for 40i, and the remaining 
kinetic energy is 4.6 lO^^^E. At the end of the prepara- 
tion, friction is turned on. 

2.5 Loading the sample 

The configurations obtained are submitted to an increas- 
ing external force along the vertical axis, whereas the hor- 
izontally applied forces remain unchanged. The vertical 
force increase is linear. 



fit) ^fo + at 



(6) 



To obtain a pressure p — fo/L that is approximately unity, 
we choose /o = 0.8/. The prefactor a determines the value 
of the very small force increment per timestep. We study 
failure in the quasi-static limit, meaning that the applied 
force at which the assembly fails is independent of a, in 
our case a = 1.28 10~^f/t. We carefully checked that 
failure is independent of a for lower values of a, but shifts 
to larger forces f{t) for higher values. At the beginning, 
a parabolic matching is applied to obtain a continuous 
differentiable force curve. 

3 Description of the simulations 

3.1 Stress-strain curve 

Fig. |2] shows a typical stress-strain curve for one simula- 
tion. Here, the deviatoric strain is defined as 



Lx L 



xO 



Ijy ±J 



2/0 



L 



xO 



L. 



yO 



(7) 



Where L^ , Ly are the horizontal and vertical length of 
the system in Fig. fll and L^o, LyQ are their initial values. 
Similarly, the deviatoric strain is 



m 






(8) 



We will refer to the slope of the curve in Fig. [2] as the 
" stiffness of the assembly" . At the beginning, the assem- 
bly is very stiff (the stiffness is comparable to /c„ = 1600, 
which is the stiffness of one contact). As the stress in- 
creases, the assembly becomes softer. The thin arrows 
indicate two "precursors" or small rearrangement events 
preceding failure. These events leave no clear sign on the 
stress-strain graph, but appear clearly when other quanti- 
ties are plotted. We will examine the precursors in detail 
in Sec. |4j Finally, the slope of the curve becomes nearly 
flat, and the assembly is very weak. Then little horizon- 
tal jumps appear that are associated with rearrangement 
processes. These are events where the walls move rapidly. 
The heavy arrow indicates the beginning of failure, which 
we define as the event where the strain crosses a threshold 
of 5%. This value is more than twice the strain just before 
failure, but much lower than the strain after failure. 

We performed several tests on systems with different 
numbers N of particles, but all having the same mass 
and size. We find that the mean value of /faii varies 
only slightly with N, whereas its variance decreases sig- 
nificantly for larger systems. This indicates that /fau is 
independent of system size. 



failure 




Figure 2: Deviatoric stress versus deviatoric strain for a 
typical quasi-static simulation with a large number of par- 
ticles (parameters defined in Secs |2.2[|2.5| and discussed in 
Sees 2.5[ 3.1). The two slim arrows indicate two precur- 
sors which are also indicated in Figs [3j [6j [7] and [8j The 
heavy arrow indicates the beginning of the failure. When 
failure happens, the strain jumps to a value close to 22%. 




Figure 3: Kinetic energy for one simulation at a — 
1.28 lQ~'^E/t. At the failure, £kin rises by many orders of 
magnitude (heavy arrow). 



3.2 Kinetic energy and vibrations 

During the simulations, the kinetic energy i?kin rises as 
the packing becomes softer due to some contacts becom- 
ing sliding and some contacts disappearing. This behavior 
is shown in Fig. [3] For instance, a two-decade increase of 
i?kin is accompanied by a one-decade decrease of stiffness. 
When the packing becomes very soft, vibrations become 
much larger in amplitude. The two precursors of from 
Fig. n\ are indicated by arrows but not visible in the fig- 
ure, for their duration is comparable to the distance be- 
tween two consecutive data points. On the other hand 
they provoke slowly damped vibrations visible in Fig. [3] 
The typical oscillation period is T = 2 lQ~^i. Finally fail- 
ure shows up as a pronounced maximum of kinetic energy 
at about t — 74t (^ 410'* oscillation periods). 



3.3 Change of volume, power injected and 
the number of contacts 

One expects the volume of the packing to decrease with 
the increase of /^xt (Fig. |4| . This decrease is very small 
for the stiff particles system (AVmax/Vb = —0.05% at e = 
0.25%). With the decrease of V, one expects the number 
of contacts M to rise. However, we observe a decrease in 
M during the same strain interval. The decrease of M 
means that the particles try somehow to avoid each other. 
It should be kept in mind that the systems we prepare are 
close to densest packing, as no friction is applied during 
preparation. Turning on friction after preparation might 
provide the reason for the decrease of M, as some particle 
motions will be blocked. After e = 0.25% the volume 
starts to increase again, while M still decreases, but at a 
lower rate. Note that the increase in V is almost linear in 
e. The rate of the decrease in M slows further down with 
e until vanishing at the failure. 

Therefore we have two regimes in Fig. HI a decrease in 
volume until e — 0.25%, and thereafter an increase in V 
until the failure. An increase in volume usually means that 
the system is loosing energy, as power has to be injected 
to decrease the volume. However, when we calculate the 
power injected at the walls (Fig. [5]) we see that is is always 
positive: the energy of the system always rises. Note that 
this behavior is unexpected: the system should simply 
explode. The key to resolving this paradox is that the 
stress is strongly anisotropic: / ~ 2/o. We will show 
that there are many qualitative differences in the systems 
behavior between the two regions in Sec. |3.5[ 

A loss in contacts M also always implies a loss of ad- 
ditional possibilities to stabilize the packing: when fewer 
contacts are present, it is more difficult for the system to 
balance a higher load. Figure [2] shows that the higher the 
actual load, the more the system deforms with a further 
increase. This evidences that a lower number of contacts 
correlates with a higher deformation. 

Besides the number of contacts, sliding contacts require 
special care. These are detailed in the next section. 

3.4 Shding contacts 

3.4.1 Number of sliding contacts 

The number of sliding contacts rises as the deviatoric part 
of the external stress on the packing is increased (Fig. [6]) . 
Previous works [HI [191 ESI [25] lead one to suppose that 
the number of sliding contacts finally attain a maximum 
at failure. However, this is not true, as the maximum 
is attained well before failure (Fig. pi). This maximum 
coincides with the minimum in volume, shown in the inlay 
of Fig|4] In small packings the existence of a maximum 
before the failure is not observed, probably because of the 
limited number of contacts. This small number allows only 
for few contact status changes at a time, and every change 
results in large changes in the stability of the granular 
assembly. Therefore very few changes can already lead to 
instability and cause failure pj]. 

We observe that when the increase in Mg becomes 
slower, the number of sliding contacts plunges at certain 
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Figure 4: Relative change of the volume (left scale) and 
relative change of the number of contacts (right scale) as 
a function of deviatoric strain. The inlay shows the same 
quantities as a function of external force. 
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Figure 6: The number of sliding contacts Mg. There 
is always a decrease of Mg before failure, followed by an 
increase at the failure. Right after the failure Mg is close 
to zero as only the contacts with the walls remain sliding 

(/^wall = 0). 



positions and then quickly recovers. The sudden plunges 
in the number of sliding contacts are the signature of pre- 
cursors. In much faster simulations (a » 10~^ in Eq. [6|, 
they are not observable. We carefully checked this on the 
basis of an equal number of data points for many simula- 
tion speeds. 

During failure, large rearrangements occur, and almost 
all sliding contacts between particles disappear or close. 
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Figure 5: Power injected into the system (left scale) and 
number of contacts (right scale) as a function of exter- 
nal force. The power injected per unit of time increases 
quickly after ///g = 1.8. Simultaneously, the number of 
contacts decreases linearly (the dashed line is a guide to 
the eye). 



3.4.2 Strength of Sliding Contacts 

Fig. It] shows the average normal force (F„) at a contact 
and compares its value to the average normal force at 
a sliding contact (i^„)siiding- (-F„)siiding is always much 
smaller than (F„), in accord with the findings of Ref. j26| : 
most sliding contacts transmit less than average forces. 
Both (Fn) and (i^n) sliding iucrcasc about linearly in the 
first half of the simulation. However, the average force 
at sliding contacts increases much faster than the aver- 
age force for all contacts. Then, later in the simulation, 
(Fn) increases faster than before, whereas the increase in 
(-Fn)siiding slows Until a plateau is reached close to the fail- 
ure. 

Note that in the beginning of the simulation the aver- 
age force at sliding contacts is very small. Therefore the 
first appearing sliding contacts are very weak contacts. 
As the deviatoric stress is increased, new sliding contacts 
appear that carry larger forces. When the maximum in 
Ms is reached, the increase in (i^jv) sliding becomes slower. 
Finally (i^„)siiding attains a maximum when Mg starts to 
decrease again. This leads to the conclusion that while 
sliding contacts constantly disappear, this does not change 
the average force at sliding contacts. This finding implies 
that sliding contacts are not part of the force chains carry- 
ing the external load. On the other hand the global mean 
(Fn) increases strongly close to the failure. We will see 
in the next subsection that this behavior is linked to the 
disappearance of contacts, so every remaining contact has 
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Figure 7: Average normal force transmitted at a contact, 
and average normal force transmitted at sliding contacts 



to carry a higher fraction of the load, on average. 



3.4.3 Understanding the evolution of Mg 

We will now examine in detail the variation in Mg. We 
will see that both the increase and the following decrease 
can be understood in terms of the frequency of transitions 
between the three contact statuses that have been intro- 
duced in Sec. |2.1| closed, sliding, and open. 

In Fig. [8] we show the most important contact status 
changes. In the beginning the most frequent transitions 
are from closed to sliding 'C— >-S'. But the number of in- 
verse transitions 'S^'C increases exponentially, and be- 
comes approximately equal to the frequency 'C— >S' at 
about ///o — 1-7, i.e., near the maximum of Ms- From 
then on, the two transitions cancel each other, i.e. their 
difference fluctuates around zero. The third most impor- 
tant transition is 'S— >-0'. Once 'S— >-C' and 'C— ^S' cancel 
each other, 'S— )-0' leads to the decrease in Mg that is 
observed in Fig. [6] before the failure. The other possible 
contact status changes (not shown in the figure) are much 
less frequent. At the failure itself at ///o ~ 2.1, all sorts of 
contact changes become very important, even those that 
are not displayed. The signature of two precursors can also 
be seen in the figure: sharp peaks in 'C— >S' and 'S— ?>C' at 
f/fo ~ 1-75 and 1.9 (See the next section for details on 
precursors) . 

Another important question is the spatial organization 
of sliding contacts, which we will investigate in the follow- 
ing. 



3.4.4 Spatial distribution of sliding contacts 
dering effects close to the failure 



or- 



To assess the spatial organization of sliding contacts, we 
investigate whether a Poisson process [27] could generate 
their observed spatial distribution. Recall that a Poisson 
process is one where a fixed number of points in a region 
are selected, with each point having an equal probability 
of being chosen. Furthermore, each point is chosen in- 
dependently of the others: the choice of point A has no 



Figure 8: Contact status changes involving sliding con- 
tacts ('S'). The heavy arrow indicates the failure. See 
Sec. 13.4.31 for details. 



influence on the probability of choosing point B. If the re- 
gion is subdivided into boxes of equal size, the probability 
of observing x points in a box is 



P{x;X) = — , x = 0,l,2, 



(9) 



Here, A is the average number of points expected in one 
box. Note that the variance of the Poisson distribution is 
equal to its mean. We will make use of this later. 

To check if the sliding contacts are distributed according 
to a Poisson process, we divide the packing into equally- 
sized square boxes of length I — Ad (d is the average parti- 
cle diameter) , and count the number of sliding contacts in 
each box. One can then compare the observed frequencies 
with the prediction in Eq. ([9|. One convenient way to do 
this is the so-called "t-test" . 

The t-test compares the variance a^ of the distribution 
to the mean Mg.box P7j : 
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(10) 



Negative i- values indicate low variance or evenness of the 
distribution while positive t-values indicate a departure 
in the direction of high variance or clumping (clustering) . 
The values in Fig. |9] show that in the beginning of the 
simulation the sliding contacts tend to be distributed ran- 
domly over the packing (t « 0). However, the very initial 
values might not be significative as the number of slid- 
ing contacts is small. As more sliding contacts appear, 
the t-value decreases and becomes negative, indicating a 
sharply peaked distribution (low variance), meaning that 
the sliding contacts repel each other: the presence of one 
sliding contact reduces the probability that a neighboring 
contact will become sliding. Later, close to the failure, 
the i-values become positive and increase strongly. This 
indicates that at the failure the sliding contacts strongly 
tend to cluster. 



Fig. 10 shows the positions of the sliding contacts just 
before failure (at f/fo = 2.15, failure is at f/fo = 2.18). 
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Figure 9: Values from the t-test (Eq.[lO|). This test checks 
the spatial distribution of sliding contacts: t < : contacts 
are more uniform than random, t > : contacts cluster. 



Figure 10: Positions of the sliding contacts for ///o = 
2.15, i.e., just before failure. At this time, the t-test yields 
i = 50. 
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Table 1: Time behavior of characteristic quantities during 
the two periods. The arrows indicate the evolution with 
increasing load (rising /^ or falling \). 



One discerns a diffuse diagonal band that crosses the sam- 
ple from lower left to upper right. At failure, a shear band 
forms in this region. The other sliding contacts are con- 
centrated at some distance in another band that is parallel 
to the shear band. The formation of the shear band will 
finally lead to the failure of the assembly. 

3.5 Two qualitatively different regimes 



As discussed in Sec. |3.3[ we can identify two regimes of 
loading with different behavior. Table [l] shows the behav- 
ior during the two periods. Most quantities in the table 
depend on the contact status changes. In the first period 
{f/fo < 1.8) the most frequent contact transition is from 
closed to sliding. Therefore the number of sliding contacts 
Ms increases with the load, while their spatial distribution 
is quite uniform (t < 0). Also the volume decreases during 
this period. In the second regime, the dominant contact 
status transition is from sliding to open (disappearing of 
formerly sliding contacts). This entails a decrease in Mg, 
and a clustering tendency (i > 0). During this period, the 
average normal force transmitted at sliding contacts does 
not increase any more. Last but not least, the change in 
volume is reversed: V increases, while the power injected 
stays positive (see Fig. [s] in Sec. 3.3). 

Towards the beginning of the second regime, precursors 
of failure appear that become more frequent with increas- 
ing /. Their appearance has been outlined in Sec. 3.4.1 
and indicated in Fig. |6] They will be discussed in detail 
in the next section. 



4 Precursors 

4.1 Definition 

Prior to the collapse of the packing, several precursors 
occur where M^ plunges and then quickly recovers. We 
define the precursor to be an event where Mg plunges by 
at least 10% of its maximum value before the failure. This 
drop in Ms varies from precursor to precursor, but the 
qualitative behavior of the precursors is always the same. 
A closer inspection in this section will show that precursors 
are initiated by an instability that gives rise to a local 
increase in the kinetic energy. The decrease in the number 
of sliding contacts is then just a consequence of the release 
of potential energy. 

Figure [6] tells us that the precursors become more fre- 
quent as the failure is approached. Therefore they might 



play an important role for the appearance of failure. In the 
next section we will see what happens at one precursor. 



4.2 Examination of a Precursor 

4.2.1 Number of sliding contacts 

To better understand the precursors, we examine in de- 
tail the first precursor indicated by the first arrow in 
Figs [2) [3) [6) [8| and |9] at f/fa « 1.71. Figure [u] shows 
Ms around this precursor. At its appearance, the drop in 
Ms is very sharp, while the recovery afterwards is slower 
and represents the relaxation to a new (force) equilibrium. 
Figure [TT] also shows the kinetic energy at that precursor. 
When Ms starts to decrease, E'kin increases quickly. How- 
ever, the kinetic energy very soon decreases again. This 
happens before Ms reaches the minimum value. The max- 
imum of i?kin can vary from one precursor to another one, 
but it is always much smaller than the maximum at fail- 
ure. This is due to the limited time in which the energy 
rises [23l. 



4.2.2 Appearance of an instability 



In Fig. 12 we show another measure of the stiffness of 



the assembly. This stiffness is calculated by reducing the 
stiffness matrix k (see appendix [A|), which contains the 
stiffnesses of all the contacts, to a scalar stiffness by mul- 
tiplication with the particle velocities v: 



k — vkv/vv . 



(11) 



Note that k contains the velocities of all the particles, 
whereas the stiffness defined earlier in Sec. 13.11 concerns 
only the walls. The advantage of k is that it can detect 
localized instabilities [23]. Specifically, fc < means that 
the packing is (at least locally) unstable, whereas A: > 
indicates that it is stable. Furthermore, k is correlated 
to the stiffness defined in Sec. |3.1| l arge positive k corre- 
spond to stiff assemblies. In Fig. |l2|we see the stability of 
the assembly at the time when the precursor appears. At 
f/fo < 1.711, the stiffness is positive. Its value does not 
change significantly until Ms starts to decrease. This hap- 
pens exactly when the stiffness becomes negative, hence 
when the packing is unstable. Shortly thereafter the stiff- 
ness becomes positive again, while Ms still continues to 
decrease. The kinetic energy rises rapidly when fc < 0. 
But at f/fo — 1.7111, k suddenly jumps to a positive 
value, and i?kin starts a rapid decline. This shows that 
the Ekin is controlled by k. We showed this dependence in 
an earlier paper on failure in small packings [23J. 

After the precursor, large vibrations appear that last 
for a long time compared to the intervals shown in Figs 
[T2l[TT1 anc(T5| These vibrations can be observed in Fig. |3] 
We anticipated in |23] that vibrations will become impor- 
tant in large packings around failure. The vibrations trig- 
gered by the precursor studied here do not cause failure, 
but Fig. [3] shows that vibrations grow as failure is ap- 
proached. These vibrations may play an essential role in 
causing the collapse of the assembly. Therefore failure 
might finally be initiated by the vibration generated by 
the precursors immediately preceding failure. 
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Figure 11: Kinetic energy at the first precursors indicated 
in Figs [2] [3j [6j |8] and [9J The energy increases quickly, and 
then decreases again. 
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Figure 12: Stability of the assembly at the occurrence of 
the precursor. When the instability appears, i?kin rises 
exponentially. The dotted horizontal line separates the 
two regions 'stability' and 'instability'. The two vertical 



lines indicate the positions of Figs. 13 14 
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Figure 13: Velocities in the assembly near the peak in Skii 
(the last point of instability at ///o = 1.71108 in Fig. 12 1 



The black particles (0.63% of the total) have above average 
kinetic energy and carry 85% of the total energy. 



Figure 14: Velocities in the assembly at ///o = 1.71112, 
shortly after the recovery of stability in Fig. [12] The black 
particles (16% of the total) have above average kinetic 
energy and carry 64% of the total energy. 



4.2.3 Localization of the kinetic energy 

When i?kin starts to rise at the beginning of the precur- 
sor, the velocities in a small region rise and become signif- 



icantly larger than everywhere else. Figure 13 shows the 
grains of the packing that carry most of the kinetic energy 
near the peak of -Ekin- But the instability for the precursor 
examined in Sec. |4.2| lasts only for a short while, and the 
velocities in this region decrease again very soon. How- 
ever, there is a wave of large movements spreading from 



this small region across the packing. Figure 14 shows the 
grains with large kinetic energy shortly thereafter, when 
the energy is propagated through the system. In this en- 
ergy spreading many different orientations are involved, 
and the propagated waves will move across the entire pack- 
ing. Looking at Figs[l2 13 and Fig. 14 we see that the 
width of the spike in Ekin in Fig. [12] is much smaller than 
the time it takes for the disturbance to cross the sam- 
ple. That means that the "high" energies (> 10~^°) occur 
only in a very localized region. This in turn shows that 
precursors are indeed "localized failure" . 



4.2.4 Where does the number of sliding contacts 
decrease? 

Figure [15] shows Ms and the i-test at the precursor. Before 
the precursor, the values are negative and do not change 
with increasing external force. When the precursor ap- 
pears, the i-values increase strongly and become positive. 
Note that the maximum positive value is much larger than 
the negative value before the precursor. It is reached at 
the time of minimum Mg. Fig. 16 shows the spatial dis- 
tribution of the sliding contacts at the maximum i-value. 



The sliding contacts disappear in some regions around the 
precursor, while the number of sliding contacts looks much 
more uniform far away from the precursor. After the pre- 
cursor, the i-values decrease again and become close to 
the values before the precursor. This indicates that order- 
ing effects appear at the precursor, and disappear again 
after the precursor 



14 and Fig. 16 



we 



Comparing Fig 
see that sliding contacts disappear in regions of elevated 
kinetic energy. Thus we conclude that the drop and subse- 
quent recovery of M^ are due to wave radiating outwards 
from the local failure. After the wave passes, the sliding 
contacts reappear, explaining why both Mg and the t-test 
return to their initial values. 

We anticipated in our last paper [23] that these local 
failures will occur. One of our findings in [23] was that 
the number of sliding contacts vanishes at the failure. This 
statement is now extended to precursors, therefore these 
can indeed be called local failures. 

4.2.5 Change in the number of contacts 

Another feature of the precursor event is a transitional 
change in the total number of contacts. Fig. [iT] shows 
the total number of contacts M as well as the number 
of sliding contacts Mg during the precursor. Surprisingly, 
these quantities are anti-correlated: The sudden drop in 
the number of sliding contacts coincides with a peak in the 
number of contacts. The contacts that are created are con- 
centrated in the high kinetic energy regions of Figs [T3| [Mj 
Therefore they are probably another effect of the wave. 
When Ms rises again, the number of contacts reduces 
slowly and attains values lower than those before the pre- 
cursor. 




1.714 



Figure 15: Values from the t-test (Eq. 10 ) and Mg during a 
precursor. The decrease in Mg corresponds to a clustering 
of sliding contacts {t > 0). 



Figure 16: Distribution of the number of sliding contacts 



in the packing at the time of the highest i- value in Fig. 15 




31120 
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Figure 17: Change in the number of contacts at one pre- 
cursor. 



Another issue raised by Fig. [17] are the permanent 
changes induced in the contact network by the precur- 
sor. For example, M decreases by about 30 between 
///o = 1.709 and 1.713 - is the precursor responsible for 
this change? Plotting M and Ms over a long time (Fig. 19 1 



shows that these changes are just part of a long term gen- 
eral trend. Furthermore, if one identifies the contacts that 
have disappeared, one finds that they are not concentrated 
anywhere in particular. 

We conclude therefore that the precursor does not lead 
to a significant change in the number of contacts. Fur- 
thermore the geometric structure of the contact network 
remains almost unchanged. Temporary changes mainly 
originate from an increase of the normal force Fpf at some 
contacts. This increase is generated by the compression 
wave radiated outward from the precursor. 

4.3 Can precursors be observed experi- 
mentally? 

While the stiffness of the assembly becomes negative and 
the kinetic energy rises, the stress-strain curve does not 
show a maximum at that time but rather a dip (Fig. |18[ ) . 
Thus in the stress-strain curve in Fig. [2] it is hard to iden- 
tify the precursor. Therefore it seems to be hard to even 
notice precursors in experimental biaxial (or triaxial) test, 
as it is difficult to detect small fluctuations of the stress- 
strain curve. 

One possibility to detect precursors in experimental in- 
vestigations is therefore to monitor the kinetic energy 
by detecting sound emissions from these regions. These 
sound emissions arise at the local grain displacements [28] . 
Sound waves of high frequency are quickly diffused [29] . 
while low frequency waves can travel the packing almost 
unchanged and can then be detected at the boundaries of 
the packing. By measuring the travel distance to differ- 
ent detectors, the spatial origin of the sound waves can be 
reconstructed [301 . However, one must take into account 
the dependence of the speed of sound on both the surface 
structure of the grains and the dimensionality (2D or 3D) 
|31j . Sound waves have been observed in triaxial tests but 
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Figure 18: Stress-strain relation during a precursor. The 



range of the data corresponds to the range in Fig. 11 At 
the precursor the relation is not linear any more, but there 
is no visible peak at that time. 



not analyzed so much. 

4.4 Summary 

The definition of precursors is based on the observation 
that the number of sliding contacts decreases suddenly at 
certain values of the external force /. A closer inspection 
of one precursor shows that there is a multifaceted behav- 
ior at this time: first of all, the precursor is initiated by 
an instability. Through this instability, the kinetic energy 
rises quickly, but only in a limited region of the packing. 
The rising energy initiates the observed decrease in the 
number of sliding contacts. The further decrease in Ms is 
then governed by wave radiation. This wave also increases 
the total number of contacts. Both of these changes are 
transitional. When the wave is gone, M and Mg almost 
return to their values before the precursor. A closer in- 
spection of the contact status changes furthermore reveals 
that only a few contacts permanently disappear (38) or 
are created (14). 

We therefore conclude that the precursors are localized 
instabilities, in contrast to failure, i.e., the global loss of 
stability. 

5 Conclusion 

We investigated the behavior of granular packings in two 
dimensions that are submitted to an increasing external 
force along one direction. We performed this biaxial test 
under quasi-static conditions where the forcing is slowly 
increased to a critical force. This critical force leads to a 
large deformation of the packing. 

We found that the time leading up to failure can be 
divided into two roughly equal periods. Many variables 
show a qualitatively different behavior in these two peri- 
ods. For example, the volume decreases in the first pe- 
riod, but increases in the second. The number of sliding 
contacts increases linearly in the first period, and then 
decreases in the second, attaining a maximum near the 



transition between the two behaviors. Thus the number 
of sliding contacts is not related univocally to the stiffness 
or the stability of the packing. The spatial organization of 
sliding contacts also changes during the simulation. Be- 
fore the maximum, sliding contacts are more uniformly 
distributed in the packing than afterwards. Comparison 
with a Poisson process shows that sliding contacts initially 
repel each other: the presence of a sliding contact reduces 
the probability that a neighboring contact will become 
sliding. After the maximum, in the second regime, the 
situation is reversed: sliding contacts are concentrated in 
specific regions. Near failure, the formation of the shear 
band is foreshadowed by a concentration of sliding con- 
tacts. These contacts cluster preferably to form diagonal 
bands. This suggests that the localization of deformation 
begins long before any shear band is visible. 

The changes in the number of sliding contacts are caused 
by the frequency of the different contact statuses transi- 
tions; in the first period, the main transition is from closed 
to sliding, leading to an increase in Ms- In the second 
period, the transition sliding to open is dominating, de- 
creasing Mg. 

Around the time of the maximum number of sliding con- 
tacts, near the transition between the two regimes, precur- 
sors begin to appear, becoming more and more frequent 
as failure is approached. These precursors are triggered 
by instabilities that lead to a sudden rearrangement of a 
small, localized number of grains who carry most of the 
kinetic energy. Precursors involve also a strong decrease in 
the number of sliding contacts, and a temporary increase 
in the number of contacts. When stability is recovered, 
the packing relaxes to a new equilibrium with properties 
(M, Ms, Skin, ■ • ■) close to the values before the precursor. 

Precursors are initiated by instabilities, therefore iner- 
tia effects become important. Hence, when investigating 
the micro-macro transition, inertia cannot be neglected 
any more, complicating the establishment of a macroscopic 
theory based on microscopic, static quantities such as the 
fabric tensor or the number of sliding contacts. 

The appearing instability during the precursor leads to 
large vibrations involving motions of all particles. These 
motions let the packing explore a larger part of the phase 
space. When approaching the critical external force, the 
packing becomes very soft, and the vibrations, triggered 
through the precursors, become larger. We therefore ar- 
gue that the precursors observed in this study should be 
significantly involved in failure. This supposition might 
be investigated in a future article. 

Financial support of the DFG through SFB716, project 
B3, is acknowledged. 

A The stiffness matrix 



In Sec. 4.2.2 we used the quantity k — vkv/vv to es- 
timate the stiffness of the packing. This stiffness is the 
sum of a mechanical part fcmoch and a geometric contri- 
bution fcgcQ. Usually the mechanical part is dominating. 



'^incch '^ "-J 



gco- 



In Ref. 



it has been shown that this 



quantity is connected to the macroscopic stiffness of the 
packing, i.e., it shows how large is the deformation for a 
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Figure 19: Number of contacts in the assembly during 
a period where three precursors appear. The loss in the 
number of contacts seems not to be connected to the ap- 
pearance of precursors. 



certain change in load (Note that in |23] k is defined to be 

"'niGchj • 

The stiffness matrix has been extensively discussed in 
Refs. 123 IS2] ■ It arises when one writes the equations of 
motion for all the particles in vector form. More specif- 
ically, one forms the vector v containing the velocities 
(translational and angular) of all the particles. In a two- 
dimensional system, it has 3A^ components, where N is 
the number of particles. If the motion is quasi-static, then 
— kv is the temporal derivative of the contact forces ex- 
erted on each particle. When the packing is stable, these 
forces balance the applied load. 

One criteria for stability is vkv > 0, i.e. positive stiff- 
ness. (We normalize vkv by vv so that fluctuations in 
velocity do not affect k. Note that — vkv/vv is the stiff- 
ness' contribution to the second derivative of the kinetic 
energy i^kin-) In small systems, failure often occurs when 
a contact status change leads to a modification of k that 
makes A: < [23] . Now vkv < implies that (the symmet- 
ric part of) k has at least one negative eigenvalue. The 
amplitude of its eigenvector grows exponentially during 
the instability. 

In large packings, the rise in E'kin is localized, there- 
fore only a small number of velocities control the change 
of k. More specifically, at least one of the eigenvalues of 
(the symmetric part of) k must be negative at the precur- 
sor, and its eigenvector v* defines which particle velocities 
grow quickly with time. This growing velocities conse- 
quently define the stiffness {k « v*kv*/v*v*). This also 
shows why the jump in k from negative to positive must 
appear exactly at the maximum in -Bkin^ at that time, the 
negative eigenvalue becomes positive. Therefore k jumps, 
as the corresponding eigenvector v* is large. Note that if 
v* did not grow so much, it would not control k. There- 
fore it is probable that small precursors with a small loss 
of sliding contacts do not lead to a negative value of k. 

Remark: the radiation of the wave is a dynamic process, 
and therefore not captured by dfc^t/dt = kv. 
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